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Dynamical models implemented on the large scale architecture of the human brain may shed 
light on how function arises from the underlying structure. This is the case notably for 
simple abstract models, such as the Ising model. We compare the spin correlations of the 
Ising model and the empirical functional brain correlations, both at the single link level and at 
the modular level, and show that their match increases at the modular level in anesthesia, in 
line with recent results and theories. Moreover, we show that at the peak of the specific heat 
(the critical state) the spin correlations are minimally shaped by the underlying structural 
network, explaining how the best match between structure and function is obtained at the 
onset of criticality, as previously observed. These findings confirm that brain dynamics under 
anesthesia shows a departure from criticality and could open the way to novel perspectives 
when the conserved magnetization is interpreted in terms of an homeostatic principle imposed 
to neural activity. 
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It has been shown that a wide class of mod¬ 
els, spanning a wide range from abstract to bio¬ 
logically detailed, reproduce large scale collective 
dynamics in the brain when they are in a critical 
regime. Here we focus on possibly the simplest 
one, the Ising model, implemented on the struc¬ 
tural architecture of the brain, and look at what 
happens when we introduce a further conserva¬ 
tion constraint: the total magnetization remains 
constant at each step. We show that this leads to 
an improved correspondence between structure 
and function at the level of modules. This phe¬ 
nomenon is increased in particular under loss of 
consciousness, when brain dynamics moves away 
from the critical regime, thus providing insights 
on how structure and function interact in the 
brain. 


I. INTRODUCTION 

One of the key challenges in the study of complex net¬ 
works is understanding the relation between structure 
and the collective dynamics stemming from it. This issue 
is of special relevance in neuroscience, where the question 
translates to how structurally distinct and distant brain 
areas dynamically interact^, both in healthy and patho¬ 
logical conditions. Recent advances in diffusion imaging 


and tractography methods allow the noninvasive in vivo 
mapping of white matter cortico-cortical projections at 
relatively high spatial resolution^, yielding a connection 
matrix of interregional structural connectivity (SC). Sim¬ 
ilarly, functional MRI can be used to obtain a functional 
connectivity (FC) matrix, by calculating the statistical 
dependencies between BOLD time series measured at dif¬ 
ferent sites of the braini^. Since the early days of con- 
nectomics, the relation between SC and FC has been a 
matter of interest, being expected but not trivial^i^. 

The intricate interplay between structure and func¬ 
tion can be investigated by simulating spontaneous 
brain activity on structural connectivity maps. Recent 
studies^^— have implemented models of dynamical os¬ 
cillators on the connectome structure^i. These compu¬ 
tational models vary from complex, biologically realis¬ 
tic spiking attractor models, describing the firing rate 
of populations of single neurons, over mean-field models 
of neuronal dynamics, down to the simple, biologically- 
nai've Ising model. All these studies agree that the 
best agreement of simulated functional connectivity with 
empirically measured functional connectivity can be re¬ 
trieved when the brain network operates at the edge of 
dynamical instability. This state corresponds to the crit¬ 
ical regime, and for the Ising model coincides with the 
maximum value of the heat capacity and of the suscepti¬ 
bility. In particular some studies showed that the resting 
activity exhibits peculiar scaling properties, resembling 
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the dynamics near the critical point of a second order 
phase transition, consistent with evidence showing that 
the brain at rest is near a critical point^^. Moreover, 
the possible origin and role of criticality in living adap¬ 
tive and evolutionary systems has recently been ascribed 
to adaptive and evolutionary functional advantages'^. 

the large-scale pattern of empirical brain correlations 
was compared with those from a large two-dimensional 
Ising model, showing that the match is optimal when 
the statistical system is close to the critical temperature. 
Remarkably, it has been recently argued that propofol- 
induced sedation and loss of consciousness move brain 
dynamics away from the critical regime^^. 

However, the Ising model on brain networks has so far 
been implemented only according to a spin dynamics in 
which magnetization is not preserved!^. Another class 
of dynamics exists, in which the total magnetization is 
preserved: it is used to describe for example alloy sys¬ 
tems, where the two different spin states naturally cor¬ 
respond to the two component atoms that comprise the 
alloyi^ and can be implemented via a pair exchange up¬ 
date rulei^. If we consider the Ising model on the human 
connectome as a model of neural activity, the conserva¬ 
tion of magnetization may be seen as a sort of homeo¬ 
static principle for the overall activity of the brain. 

The question we address here is whether an Ising model 
with conserved magnetization on the human connectome 
would be a more suitable model for the functional con¬ 
nectivity of brain, in particular under anesthesia, where 
previous works have hypothesized a departure from crit¬ 
icality. 

Under anesthesia, the brain spans a dynamical reper¬ 
toire that is reduced with respect to wakefulness. This 
would result in an increased correspondence between 
structural and functional connectivity^i^. Following the 
reasoning and the results we think that this corre¬ 

spondence is to be sought at the level of modules rather 
than at the level of individual links. 


II. METHODS AND DATA 

A. Dynamical Ising models on brain networks 

The study of the Ising model, an appealing descrip¬ 
tion of phase transitions in ferromagnets, played a fun¬ 
damental role in the development of the modern theory 
of critical phenomena (see, e.g.^ and the recent review 
for Ising’s model 90*^ birthday^!). In the original model, 
a regular lattice is populated by 2-state spins, assuming 
one of the two values ai = ±1. Pairs of nearest neigh¬ 
bours spins interact so as to favour their alignment. The 
Hamiltonian of the system is given by 

H = 

<b) 

the sum being over nearest neighbours pairs on the lat¬ 
tice, the positive coupling J favouring ferromagnetic or¬ 
der. For spatial dimension d > 1 the model exhibits, in 
the thermodynamic limit, a phase transition with finite 


critical temperature Tc, such that above Tc the spatial ar¬ 
rangement of spins is disordered with an equal number of 
up and down spins. Below Tc the magnetization is non¬ 
zero and distant pairs of spins are strongly correlated. 
All the equilibrium properties of the Ising model can be 
obtained from the partition function Z = exp{—j5'iV), 
the sum being over the configurations of the system and /3 
being the inverse temperature. Dynamical rules leading 
to the same equilibrium are not unique, all the possibil¬ 
ities being fixed by the detailed balance condition. Fun¬ 
damentally, the spin dynamics may or may not conserve 
the total magnetization, depending on whether the Ising 
model is being used to describe alloy systems, where the 
magnetization is conserved as it is related to the composi¬ 
tion of the material, or spin systems where magnetization 
is not conserved. 

According to Glauber dynamics^^ the magnetization 
is not conserved and each spin is sequentially considered 
and flipped with probability PfUp = (1 + exp{/3AE))~^, 
where AE is the energy difference associated to the 
spin flip. Here we consider the Kawasaki spin-exchange 
dynamicsi^ which conserves the magnetization: two 
spins randomly chosen are swapped with probability 

exp(-AF;), 

where AE is the variation of the energy corresponding 
to exchanging the two spins. A full iteration consists 
in tentatively updating all the spins (pairs of spins) for 
Glauber (Kawasaki) dynamics. 

Turning now to couplings, let us denote Atj the sym¬ 
metrical structural connectivity matrix. The Hamilto¬ 
nian of the Ising model on the network is 

H — ^ ) Jij O'iO'j , 

where the couplings are given by Jij = (3Aij and the pa¬ 
rameter P plays the role of an inverse temperature. Since 
we deal with finite size systems, they exhibit a (pseudo)- 
transition between the disordered phase and the ordered 
one, corresponding to the peak of the specific heat (and 
of the susceptibility, for the case of Glauber dynamics). 

In2^ the Ising model with Glauber dynamics was imple¬ 
mented on the human connectome matrix oSii at two dif¬ 
ferent spatial scales, 998 and 66 nodes, and the directed 
and undirected information transfer between nodes was 
then quantified. 

Spin correlations were evaluated using the classical lin¬ 
ear Pearson correlation. The pairwise transfer entropy 
TE, measuring the information flow from spin i to spin 
j in each pair connected by a link in the underlying net¬ 
work was computed as follows: 

TErj= Y. E E 

CTj—il Sj—il Sj —ihl 

where p{Tj,Ti) is the fraction of times that the config¬ 
uration (Ej,Si) is observed in the data set, and similar 
definitions hold for the other probabilities. 
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It was shown that at criticality the model displays 
the maximal amount of total information transfer among 
variables, with patterns consistent for both the coarser 
and the denser scale. Given the fact that in this line 
of research we are particularly interested in information 
flow in networks, we have verified that the peak of the 
information transfer corresponds to the peak of the spe¬ 
cific heat; therefore a range of temperatures around the 
peak of the information flow (or, equivalently, the peak 
of the specific heat) will be taken as the critical regime 
of the system. It is worth mentioning that it has been 
observed in the regular 2D lattice Ising model that, dif¬ 
ferently from the pairwise information flow which peaks 
at criticality, the global information transfer peaks on 
the disordered side of the transition, and the asymmetry 
observed iii 2 & as also been observed in^ when consid¬ 
ering both Ising and Kuramoto (from the viewpoint of 
correlations in phase synchronisation). 


B. Data 

The fMRI data that we consider in this work were 
recorded from healthy subjects in awake conditions and 
during propofol anesthesia. The motivation for the study, 
the underlying physiological issues, and the protocol are 
extensively described in^S. The functional MRI (fMRI) 
data was preprocessed with FSL (FMRIB Software Li¬ 
brary v5.0). The first 10 volumes were discarded for 
correction of the magnetic saturation effect. The re¬ 
maining volumes were corrected for motion, after which 
slice timing correction was applied to correct for tem¬ 
poral alignment. All voxels were spatially smoothed 
with a 6 mm FWHM isotropic Gaussian kernel and af¬ 
ter intensity normalization, a band pass filter was ap¬ 
plied between 0.01 and 0.08 Hz. In addition, linear and 
quadratic trends were removed. We next regressed out 
the motion time courses, the average CSF signal and 
the average white matter signal. Global signal regres¬ 
sion was not performed. Data were transformed to the 
MNI152 template, such that a given voxel had a vol¬ 
ume of 3mm X 3 mm x 3mm. Finally we obtained 
116 time series, each corresponding to an anatomical 
region of interest (ROI), by averaging the voxel sig¬ 
nals according to an anatomical template^^. We se¬ 
lected this partition for being the most used in fMRI 
connectivity analysis, and because it includes subcorti¬ 
cal structures. For the diffusion MRI data, we used the 
publicly available data contained in the Nathan Kline 
Institute- Rockland sample described and downloadable 
at http: / / fcon.lOOO. projects. nitre. org/ indi/ pro /eNKI_ 
RS_TRT/Frontpage.html. As a first step, the images 
were corrected for motion and eddy currents due to 
changes in the gradient field directions of the MR scan¬ 
ner. In particular, the eddy-correct tool from FSL was 
used to correct both artifacts, using an affine registra¬ 
tion to a reference volume. After this, DTIFIT was used 
to perform the fitting of the diffusion tensor model for 
each voxel. Then, a deterministic tractography algorithm 
(FAGT)^ was applied using TrackVis^, an interactive 
software for fiber tracking. Two computations were per¬ 


formed to transform the anatomical atlas to each individ¬ 
ual space: (I) the transformation from the MNI template 
to the subject’s structural image (TI), and (2) the trans¬ 
formation from the subject’s TI to the subject’s diffusion 
image space. Combining both transformations, each at¬ 
las region is transformed to the subject’s diffusion space, 
allowing to count the number of reconstructed stream¬ 
lines connecting all ROI pairs. 

It is worth to note that for group level analyses at the 
scales considered in this study it is not relevant that the 
structural connectivity used for the simulations is not 
the one obtained from the same subjects for which the 
functional connectivity was computed and compared^. 


C. Modularity and modular similarity between networks 

A key concept in network theory is modularity^. It 
describes how efficiently a network can be partitioned 
in sub-networks, and it is particularly relevant when it 
comes to study the interplay between anatomical segre¬ 
gation and functional integration in the brain^. Maxi¬ 
mization of the modularity Q allows to identify commu¬ 
nities. Here we use the algorithm described in^^ which 
also takes into account negative weights, a situation that 
frequently arises in functional networks. The resolution 
parameter is set to its default value 7 = 1 and, for 
each network, the algorithm was run 1000 times choos¬ 
ing the maximal Q and the corresponding partition. In 
order to compare, at the modular level, two networks 
with the same nodes, we calculate the similarity of their 
partitions quantifying it by the mutual information ap¬ 
proach described in^. The code used to compute these 
quantities is contained in the Brain Connectivity Toolbox 
(https://sites.google.com/site/bctnet/) 


III. RESULTS 

We implemented the Ising model on the structural 
brain networks, and evaluated the transfer entropy and 
the spin correlations, as described in the previous sec¬ 
tion. The results shown here correspond to the average 
over 1000 runs, each consisting of 30000 full iterations of 
the lattice (after discarding the transient). Since we deal 
with a small system, and we are not interested in the low 
temperature limit, we always assumed zero magnetiza¬ 
tion for Kawasaki dynamics, i.e. the starting configura¬ 
tion consisted in an equal number of plus and minus spins 
randomly assigned to nodes. In other words, due to the 
small size, we assume vanishing equilibrium global mag¬ 
netization. We start considering, in figure o, the follow¬ 
ing problem: to which extent are the functional patterns 
of the Ising systems shaped by the underlying topology, 
at the level of individual links? We compare the anatomi¬ 
cal network with the functional networks provided by the 
dynamical system, as a function of the inverse tempera¬ 
ture. The following quantities are depicted, derived from 
the Ising model with Kawasaki dynamics on the struc¬ 
tural connectome corresponding to awake conditions: the 
Pearson correlation between the transfer entropy TE and 
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the coupling J over all the anatomically connected pairs, 
and the Pearson correlation between spin correlations c 
and J over all the anatomically connected pairs. The 
plots of the average absolute value of the correlation and 
the average transfer entropy, both peaking at criticality, 
are also displayed for comparison. As figure m shows, 
in the temperature regime at which information trans¬ 
fer and correlations are higher (identified with the crit¬ 
ical regime, see below), the link correlation of TE and 
c, with J, is minimal; in other words the critical states 
appear to be the ones at which the functional pattern is 
minimally shaped by the details of the underlying struc¬ 
tural network. These findings are in line with a previous 
study on Ising models implemented on the connectome, 
in which the best fit between model and empirical cor¬ 
relations was observed when the entropy of the model 
attractors started to increase, and not at its peabi^. 
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FIG. 2. Top: matrices corresponding to empirical data. 
Structural connections (left), average correlations in wake¬ 
fulness (center), average correlations in anesthesia (right). 
Bottom: average correlations in simulated time series for 
three values of the inverse temperatures, corresponding to 
relevant points in the curves of the following figures. 

As depicted in figure (jS)), when the link-wise correla¬ 
tion between model and empirical functional patterns is 
considered, the match between model and empirical cor¬ 
relations is higher in wakefulness for Glauber dynamics, 
and in anesthesia for Kawasaki dynamics, in the respec¬ 
tive critical regimes for each dynamics. Moreover, under 
sedation, the Kawasaki dynamics results in a better fit 
and in a clearer separation between wake and anesthesia, 
compared to the Glauber dynamics. 


FIG. 1. The Ising model with Kawasaki dynamics is imple¬ 
mented on the 116 regions structural connectome. The fol¬ 
lowing quantities are depicted as a function of the inverse 
temperature /3: the correlation between the value of TE and 
J over all the anatomically connected pairs (red line, with 
bullets); the correlation between the value of c and J over all 
the anatomically connected pairs (cyan, dotted line); the nor¬ 
malized average absolute value of the correlation (blue, full 
line); the normalized average transfer entropy (green, dashed 
line). 

Next, in order to elucidate whether and how the fit be¬ 
tween model and empirical correlation change with the 
level of consciousness, we consider both the fMRI data 
recorded from healthy subjects in awake conditions and 
during propofol anesthesia, as well as the correspond¬ 
ing structural data. Varying the temperature, we im¬ 
plemented the Ising model with Kawasaki dynamics on 
the structural architecture connecting the 116 ROIs. We 
compare the corresponding spin correlations with the em¬ 
pirical functional correlations. 

In order to visualize the typical patterns of the connec- 
tomes under exam, in figure ([2]) we depict the 116 x 116 
structural connectivity matrix, the empirical functional 
connectivity matrix for wake and anesthesia conditions, 
and the patterns of correlations of the Kawasaki model 
tuned at three different temperatures corresponding to 
relevant regimes in the curves described below (greatest 
linkwise correlation, greatest mean squared error, max¬ 
imum mutual information between structural and func- 


Glauber Kawasaki 





FIG. 3. The link-wise correlation between the model spin 
correlations and the empirical functional connectivities is de¬ 
picted as a function of the inverse temperature /3 for wake¬ 
fulness (dashed red line) and anesthesia (full blue line). Left 
panel: Glauber dynamics. Right panel: Kawasaki dynamics. 

In figure (jH), the match between empirical and spin 
correlations is measured in terms of the mean square 
error between the two patterns, i.e. the average of 

( spin empirical\‘^ n . r v ■ • spin 

CjJ — Cj^- 1 over all pairs ot brain regions, cT 

being the spin correlation of the Ising models and 
^empirical empirical functional connectivity. Results 
show that again, the best match for anesthesia is bet¬ 
ter than that for wake conditions, both for Glauber and 
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Kawasaki dynamics. Using the mean square error, how¬ 
ever, Glauber dynamics showed a better fit. 


Glauber 



Kawasaki 



/3 
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FIG. 4. The mean square error between the model spin corre¬ 
lations and the empirical functional connectivities is depicted 
as a function of the inverse temperature 0 for wakefulness 
(dashed red line) and anesthesia (full blue line). Left panel: 
Glauber dynamics. Right panel: Kawasaki dynamics. 
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FIG. 5. The mutual information between the partitions of 
the empirical and model functional networks, as a function 
of the inverse temperature 0 is depicted for Kawasaki and 
Glauber dynamics, and for wake and anesthesia conditions. 
The modular decomposition is obtained by maximization of 
modularity. 


IV. DISCUSSION AND CONCLUSIONS 


However, it has been shown that a modular compari¬ 
son is better suited to investigate the interplay between 
structure and function^i. Hence, in order to better de¬ 
scribe the relation between empirical and simulated func¬ 
tional patterns at the modular level, we have evaluated 
the mutual information between partitions (obtained by 
maximizing the modularity) as a function of 0 to quan¬ 
tify the relation between empirical functional correlations 
and spin correlations from the model, for wake subjects 
and subjects under anesthesia. In figure ([5]) the mutual 
information is plotted for both Kawasaki and Glauber 
dynamics. In wakefulness, where the same structural 
connectivity subserves a wide repertoire of activity, we 
observe a reduction of the mutual information in the criti¬ 
cal regime with respect to the disordered phase, both for 
Glauber and Kawasaki dynamics. This behavior is re¬ 
versed, leading to increased mutual information between 
structural and functional models, for the Kawasaki dy¬ 
namics under sedation. 

These results speak to the fact that Ising models tuned 
at criticality result in connectivity matrices with a gener¬ 
ally good linkwise resemblance with the empirical ones, 
comparable or even better than the one obtained with 
more biologically precise models^i^^. This similarity is 
represented by a peak in the curves in figures [3] which fol¬ 
lows a trough, corresponding to the peak in mean squared 
error of|4l According to the pairwise metric though, the 
best resemblance remains the not-so-interesting one cor¬ 
responding to the limit of high temperatures. On the 
other hand, if we look at the mutual information between 
structural and functional modules, we can observe that, 
with respect to the disordered phase, Kawasaki dynamics 
in the critical regime leads to a decreased match between 
structural and functional modules in wakefulness, and an 
increased one in anesthesia. These results, in line with 
our hypotheses, are not evidenced by Glauber dynamics. 


We have considered pair exchange update rules for the 
Ising model^ implemented on the structural brain net¬ 
work at the macroscale, using a data set of healthy sub¬ 
jects scanned during quiet wakefulness and during deep 
sedation, a condition in which the structure-function re¬ 
lation is modified. Our results show that the structure- 
function relation is strengthened under anesthesia, both 
at the link and modular level, compared to wake con¬ 
ditions. Having shown that at criticality the functional 
pattern is less dependent on the underlying structural 
network, it follows that anesthesia takes the brain dy¬ 
namics farther from the critical regime, in accordance 
with previous evidence. 

Moreover, at the modular level we obtained a bet¬ 
ter match with empirical functional correlations using 
Kawasaki dynamics compared to the more common 
Glauber dynamics for which the magnetization is not 
preservediS. This improved match suggests that the 
conservation law of the Kawasaki dynamics might ad¬ 
mit a physiological counterpart. A possible interpreta¬ 
tion is seeing it as an effective implementation, due to 
time scales separation, of the coupling from metabolic re¬ 
sources to neural activity, a key ingredient that is missing 
in neural models on the connectome^. 

In agreement with recent theoretical frameworks^, our 
results suggest that a wide range of temperatures corre¬ 
spond to criticality of the dynamical Ising system on the 
connectome, rather than a narrow interval centered in a 
critical state. In such critical conditions, the correlational 
pattern is minimally shaped by the underlying structural 
network. It follows that, assuming that the human brain 
operates close to a critical regime^, there is an intrin¬ 
sic limitation in the relationship between structure and 
function that can be observed in data. We have shown 
that empirical correlations among brain areas are bet¬ 
ter reproduced at the modular level using a model which 
conserves the global magnetization. The most suitable 
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way to compare functional and structural patterns is to 
contrast them at the network level, using, e.g., the mu¬ 
tual information between partitions like in the present 
work. 

During the awake resting state, spontaneous brain 
activity constantly fluctuates across brain regions, ex¬ 
hibiting a rich repertoire of functional connectivity pat¬ 
terns. Previous studies accounted for long-range resting- 
state functional connectivity persisting even after loss 
of consciousness^^— . A recent study on monkeys^^ 
posited that the role of structural connectivity in sculpt¬ 
ing functional connectivity maps changes during wake¬ 
fulness and anesthesia. According to the authors, wake¬ 
fulness seems to be characterized by a rich repertoire 
of connectivity patterns, while the functional connectiv¬ 
ity patterns under sedation follow the underlying brain 
structure. Another study reports increased similarity be¬ 
tween whole-brain anatomical and functional connectiv¬ 
ity networks during deep sleep^S. Our results, showing 
that the structure-function correspondence is enhanced 
under anesthesia at the modular level, are in accordance 
with this evidence. 

Summarizing, we have considered the Ising model as a 
valuable tool to explore large scale brain dynamics. We 
have shown that the conservation of magnetization leads 
to better correspondence between structure and function, 
notably where this is more expected, that is in loss of 
consciousness, again speaking to a less critical dynami¬ 
cal regime. The reason of this improvement might lie in 
the features added by the Kawasaki dynamics to large 
scale connectivity (i.e. negative correlations). Its bio¬ 
logical counterparts could be intuitively found in some 
homeostasis mechanism or metabolic constraint, but no 
validation tools for this conjecture exist at the moment. 

Finally, our results confirm that matches between 
structure and function should be sought at the modu¬ 
lar level rather than among individual links, in line with 
recent wor k^^i^^d^ . 
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